Fixation in Evolutionary Games under Non- Vanishing Selection 
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One of the most striking efTect of fluctuations in evolutionary game theory is the possibility for 
mutants to fixate (take over) an entire population. Here, we generalize a recent WKB-based theory 
to study fixation in evolutionary games under non- vanishing selection, and investigate the relation 
between selection intensity w and demographic (random) fiuctuations. This allows the accurate 
treatment of large fluctuations and yields the probability and mean times of fixation beyond the 
weak selection limit. The power of the theory is demonstrated on prototypical models of coopera- 
tion dilemmas with multiple absorbing states. Our predictions compare excellently with numerical 
simulations and, for finite w, significantly improve over those of the Fokker-Planck approximation. 
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INTRODUCTION 

Evolutionary game theory (EGT) provides a natural 
theoretical framework to describe the dynamics of sys- 
tems where successful types or behaviors, as those aris- 
ing in biology, ecology and economics |0, |^ , are copied by 
imitation and spread. Evolutionary stability is a crucial 
concept in EGT and specifies under which circumstances 
a population is proof against invasion from mutants 0, || . 
This notion was shown to be altered by finite-size fluctu- 
ations and led to the key concept of evolutionary stability 
in finite populations The latter is closely related to 
the notion of fixation |j, referring to the possibility 
for mutants to take over (fixate) an entire population of 
wild species individuals. Furthermore, evolutionary dy- 
namics is characterized by the interplay between random 
fluctuations [|| and selection, that underlies adaptation 
in terms of the different reproduction potential (fitness) 
of the individuals. Thus, a parameter was introduced to 
measure the selection intensity In this context, the 
fixation probability of a species has been calculated for a 
finite two-species population in the weak selection limit 
of vanishingly small selection intensity |[ |^. This 
limit is often biologically relevant and greatly simplifies 
the analysis (treating selection as a linear perturbation). 
However, the behaviors obtained under strong and weak 
selection are often qualitatively different (see e.g. [|[ D). 

In this Letter, we study fixation under non-vanishing 
selection in EGT and provide a comprehensive analy- 
sis of the combined influence of non- vanishing selection 
and random fluctuations. As exact results for the fix- 
ation probability and mean fixation times (MFTs) are 
rarely available and often unwieldy (see e.g. [|[ ^, |7|), our 
analysis relies on the WKB (Wentzel-Kramers-Brillouin) 
approximation method |^ directly applied to the un- 
derlying master equation j^. This technique was re- 
cently used to treat generic birth-death systems that un- 
dergo metastable switching or extinction Im- 
portantly, here we generalize the WKB formalism to sys- 
tems with multiple absorbing states. This theory ac- 



curately accounts for the large fluctuations not aptly 
captured by the Fokker-Planck approximation 

(FPA) 0. We illustrate our method on two classes of 
prototypical models of cooperation dilemmas, the anti- 
coordination and coordination games, where a coexis- 
tence state separates two absorbing states in which the 
population is composed of only the fixated species while 
the other goes extinct ||l|, ^. We compute the fixation 
probabilities, the MFTs, as well as the complete prob- 
ability distribution function (PDF) of population sizes, 
and show that our theory is superior to the FPA for fi- 
nite selection strength. 



THE MODELS 

In EGT, the fitness, or reproduction potential of an 
individual, is determined by the outcome, called payoff, 
of its interaction with the others as prescribed by the 
underlying game [Q. In fact, when two A— individuals 
interact, both receive a payoff a. If an individual of type 
A interacts with another of type B, the former receives 
6 while the latter gets a payoff c. Similarly, when two 
B— individuals interact, both get a payoff d. Now, as- 
sume that in a population of size N there are n individ- 
uals of type A ( "mutants" ) and — n of type B ( "wild 
type"). The respective average payoffs (per individual) 
are IVA{n) = {n/N)a + [{N-n)/N]b and IlB{n) = 
{n/N)c+ [{N — n)/N] d while the population mean 
payoff is n(n) ^ [nIlA{n) + {N ~ n)nB(n)] /N. For in- 
finite {N — > oo) and well-mixed populations, the density 
X = n/N of the A species changes according to its rela- 
tive payoff and obeys the replicator dynamics, given by 
the rate equation |^ 



xiiiA - n). 



(1) 



Here, we are particularly interested in anti-coordination 
games (ACG), where c> a and b > d, and in coordina- 
tion games (CG), where a > c and d > b. In addition to 
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the absorbing states n = and n = N, ACG and CG ad- 
mit an interior fixed point associated with the coexistence 
of A and B species at a density x* = {d — b)/{a — b — c + d) 
of A's. According to the rate equation (|^), x* is an at- 
tractor in ACG and a repellor in CG, whereas x — and 
X — 1 are repeUing fixed points in ACG and attracting 
in CG. 

To account for fluctuations arising when the popula- 
tion size is finite, the evolutionary dynamics is imple- 
mented in terms of fitness-dependent birth-death pro- 
cesses ^ describing, e.g., the evolution of the prob- 
ability Pn{t) to have n individuals of type A at time t: 



dPnjt) 
dt 



- "^n-lPn-l + Tn+i-Pn+l - \Tn + ]Pn- (2) 



Here, an individual chosen proportionally to its fitness 
produces an identical offspring which replaces a randomly 
chosen individual and the total population size N 
is conserved. Thus, in the master equation (H), the re- 
action rates for the birth/death transitions n — > n ± 1 
are given by r± = x^UA{n)jB[n))n{N ~ n)/N^, 
where x^(n) are functions of the fitness of each species, 
Ja {n) = 1 — w + wWa (n) and fs (n) — 1 ~ w + wll b {n) . 
As often in EGT, we focus on systems evolving ac- 
cording to the fitness-dependent Moran process (fMP) 
for which x^{n) = /^[(n/A^)/^ + (1 - n/N)fBr^ and 
X'{n) = fB[{n/N)fA + (1 - n/N)fB]-^ § p|. It is 
worth noticing that X~''(n) and x i^) intersect only at 
the fixed point value n ~ Nx* for < ti < A^, which en- 
sures that the properties of the replicator dynamics (|l|) 
are recovered when N ^ oo |l^ . 

The fitnesses Ja [n) and Jb (n) are comprised of a base- 
line contribution [the constant (1 — w)] and a term ac- 
counting for selection [wll^ for /a], where the parameter 
< w < 1 measures the selection intensity ^. The 
latter is weak for w — !• 0, when T± oc n{N - n)/N'^, 
and strong for it; — > 1, when the baseline fitness becomes 
negligible. As n e [0, A^] and n — 0,N are absorbing, the 
boundary conditions to Eq. (||) are = = 0. 



WKB THEORY OF ANTI-COORDINATION 
GAMES 



Our WKB-based approach is presented in the frame- 
work of ACG {e.g. snowdrift and hawk-dove games 
where the absorbing states n = or a; = (all B s), 
and n = A^ or X = 1 (all A's) are separated by the 
interior attractor x* [in the language of the rate equa- 
tion (|l])]. However, in the presence of noise x* becomes 
metastable, which is very naturally accounted by our the- 
ory. For A^a;* ^ 1, after a short relaxation time tr, the 
system settles into a long-lived metastable state whose 
population size distribution is peaked in the vicinity of 
A^a;* p3t|. This implies that fixation of either species 



occurs only in the aftermath of a long-lasting coexis- 
tence. At t ^ tr, only the first excited eigenvector of 
(^, 7r„, called the quasi-stationary distribution (QSD), 
has not decayed and hence determines the shape of the 
metastable PDF. Indeed, at t 3> the higher eigen- 
modes in the spectral expansion of Pn{t) have already 
decayed, and the metastable dynamics of the population 



sizes PDF satisfies [13 

Pn{t) 



7r„e 



t/r 



for ne[l,N-l], 



(3) 



where J2n ""n = 1- Thus, ai t ^ tr the dynamics of the 
probabilities to be absorbed at n = and n = N satisfies 

Po{t) 0(1 - e-*/-) , Pjv(t) ^ (1 - 0)(1 - e-*/^). (4) 

Here, 4>^ = (f) and (j)^ = 1 — (j) are the fixation proba- 
bilities of the B and A species, respectively, r is the un- 
conditional MET, and a very strong inequality r 3I> tr 
holds. The fixation probability and MET are determined 
by the fluxes into the absorbing states. Therefore, using 
Eqs. @) and (jj), one obtains 

T = [Tf TTi + T^_^iiN-i] ^ , and = Tf ttit. (5) 

Similarly, the respective conditional METs of species 
A and B (conditioned on the flxation of type A and 
B, respectively) are — [T^ -^tttv-i] and = 

[Tj^tti] ^. According to Eq. (||), these quantities are de- 
termined once we have obtained tti and i:n~i from the 
full expression of the QSD that we now compute. 

The QSD satisfies the quasi-stationary master equa- 
tion, obtained by substituting Eq. (^) into (|^) and ne- 
glecting the exponentially small term 7r„/T (to be verified 
a posteriori): 



„_i7r„-i 



T~ 



-1 - 



T-]7r„ = 0. (6) 



For A^ ^ 1, we define the transition rates T±{x) — 
[0 ^ continuous functions of x, and treat Eq. (|^) 
by employing the WKB ansatz [p|-|ri]| 

7r„ = 7rj;7v = 7r(x) = Aeyip[-N S{x) - Si(x)] , (7) 

where S{x) and Si{x) are respectively the system's action 
and its amplitude, and ^ is a constant prefactor intro- 
duced for convenience. The WKB approximation is here 
an asymptotic series expansion in powers of 1/N based 
on the exponential ansatz (0) (see, e.g., [||-0) (l8|. Sub- 
stituting (0) into Eq. (||) yields closed equations for S{x) 
and 5*1(3::). To leading order, similarly as in Hamiltonian 
systems, the action obeys the Hamilton- Jacobi equation 
H(x, S') — 0. In this case, the underlying Hamiltonian 



H{x,p) = T+{x){eP - 1) + T-{x){e-P - 1), 



(8) 



where we have introduced the auxiliary momentum 
p{x) = dS/dx S^l^. Therefore, to leading order. 
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the "optimal-path" followed by the stochastic system, 
from the metastable state to fixation, is Pa{x) = 
— In [T+{x)/T-{x)], corresponding to the zero-energy tra- 
jectory H{x,Pa) = with non-zero momentum. The ac- 
tion along Paix) is 



S{x) 



in[r+(e)/r-(c)]de 



(9) 




-15r' 



Performing the subleading-order calculations, one ob- 
tains Si{x) = {l/2)ln[T+{x)T-{x)] 0,0. Imposing 
the normalization of the Gaussian expansion of the QSD 
about X — X* , one finds the constant A, yielding 
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7:{x) = T+{x*) 



S"{x*) 



2tiN T+{x)T-^{x) 



'N[S{x)~S(x')] 



(10) 



This expression is valid sufficiently far from the bound- 
aries, where T±{x) = 0(1) [0, and generally leads to a 
non-Gaussian QSD with systematic deviations from the 
Gaussian approximation near the tails, as illustrated in 
Fig. g(a). 

To obtain the full QSD we need to match the WKB 
result (^0|) with the solution of Eq. (^ in the vicin- 
ity of the absorbing boundaries, where the transition 
rates can be linearized [ETl. For instance, near x — 
T±{x) ~ xr^(O), so Eq. (|) yields (n - l)r;(0)7r„_i + 
(n + l)ri(0)7r„+i - 7i[r^(0) + r_l(0)]7r„ = 0. Its recur- 
sive solution is 7r„ = (7ri/7i)(i?Q — l)/(i?o — 1), where 
i?o — Tl{Q)/TL{0). Matching this expression with the 
leading order of Eq. (fol) in the vicinity of a; = yields 



FIG. 1: (Color online), (a) InTTp vs. n (with N = 150): 
theoretical predictions [Eqs. (|lo|)-(|l3|)] (solid) compared with 
numerical results (dashed) and with the Gaussian approxi- 
mation of the QSD (dashed-dotted). (b) lnr~^ as a function 
of A^: theoretical predictions [Eqs. (|), (|l^)-(|l§)] (solid) and 
numerical results (symbols). Parameters are a — 0.1, b = 0.7, 
c = 0.6, d = 0.2, w = 0.5 and the system follows the fMP. 
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(11) 
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A similar analysis at a; ~ 1 with Ri = 7^(1)/7^(1) gives 
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/ ATQ'U^* \ n' (n~*\('D A\ FIG. 2: (Color online), (a) Inr vs w. theoretical [Eqs. 051), 

^N-l = \l J- ' ^+ i^^n^i-l )^M5(.-)-5(i)]. (12) (0)-(0)] (solid) and numerical results (symbols), (b) D^- 
zvr \/7^(l)7L(l) pendence of Inr on the initial number n of A's, for w — 0.2, 

0.5 and 0.8 (bottom to top): comparison between theoret- 
ical (solid) and numerical (dashed) results, (c) Theoretical 
[Eq. (114)] (solid) and numerical (symbols) results for the ra- 
tio 0^ vs w. (d) Same as in panel (b) for (p^/cf)^ {w 
grows from top to bottom). Parameters are a = 0.1, b = 0.7, 
c = 0.7, d = 0.2, N = 200 and the system follows the fMP. In 
the numerical results of (a) and (c), n is chosen sufficiently 
large so that fixation does not occur immediately (see text). 



Hence, the expressions (^0|)-(^2|) provide us with the com- 
plete QSD. 



FIXATION IN ANTI-COORDINATION GAMES 



We now apply the general results obtained in the pre- 
vious section to study fixation in ACG, when the system 
follows the fMP. In this case the action given by Eq. (||) , 
becomes 

S{x) = [B/{B -A)-x] \x)[Ax + B{1 - x)] 

+ [D/{C -D) + x]\x)[Cx + D{l-x)], (13) 

where A — 1— w+u>a, B = 1—w+wb, C = l—w+wc, and 
D = 1-w + wd |T|. Provided that A^[5(l) - ^(a;*)] > 1, 
and A^[5(0) — '5'(a;*)] 1 (which imposes a lower bound 



on w), the MFTs and fixation probability are obtained 
from Eqs. (|) and (|ll|)-(|l|) with Tf = T+_^ ~ A^-i. 
These results generalize those obtained previously in the 
limiting cases Nw ^ 1 || and w = 1 (for which 
A = a, B = b, C = c, a,nd D = d) ^. As illus- 
trated in Fig. |l](b), one finds that the unconditional 
MFT asymptotically exhibits an exponential dependence 
on the population size N, t (x A^i/^e^'-^"'^'^ ^\ where 
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the governing exponent S = min [S'(O), S'(l)] is readily 
obtained from (|l|). With @-(|l|), this confirms that 
7r„/r is indeed exponentially small. For < w < 1, one 
finds that S increases monotonically with w, as shown 
in Fig. ||(a). Here (as in our other figures), the theo- 
retical predictions are compared with the numerical so- 
lution of the master equation yielding an excellent 
agreement. It also follows from (H),(|ri|)-(|l^) that for 
^ 1 and small (but not too small) selection inten- 
sity, ^ w ^ 1, the conditional MFTs grow expo- 
nentially as - ]^l/2^Nw{a~cf/[2{c^a+b-d)]^ ^^^^ _ 
j^l/2^Nw(b~df/[2(c-a+b-d)]^ ^-^^Yi T = T^T^ /{t'^+T^) 

min(r"*, t^). 

As our approach assumes that fixation occurs after the 
metastable state is reached, the expressions obtained for 
the MFTs are independent of the initial number n of A's, 
when n 3> 1. This is confirmed in Fig. ^b) where theory 
and numerical results agree excellently. 

The ratio 0^/0^ = — 1 between the fixation prob- 
abilities of the A's and B's allows to understand the in- 
fluence of selection and the interplay between selection 
and demographic stochasticity. Indeed, with Eqs. (H), 
and ([ll])-([T^), our theory yields 

In Fig. 11(c), we show the ratio (j)^/(f>^ and find a nontriv- 
ial exponential dependence on w in excellent agreement 
with numerical calculations. Contrary to the neutral case 
w = (not covered by our theory) , where the ratio of fix- 
ation probabilities strongly depends on the initial number 
of A's @, Eq. (|l|) predicts that 0^/0^ is independent 
of the initial condition when the selection strength is fi- 
nite. Indeed, the numerical results presented in Fig. ^(d) 
confirm that for n ^ 1, the ratio (f)^ / (p^ coincides with 
( p^ and becomes independent of n when w is nonzero 
(for w the convergence requires n ~ Nx*). 

WKB THEORY AND FIXATION IN 
COORDINATION GAMES 

As a further illustration of our theory, we accurately 
compute the fixation probability in CG {e.g. stag- hunt 
game [|l|). Here, the fixed point x* is a repellor while 
X = 0, 1 are attracting, hence there is no metastability 
and fixation occurs quickly . As a result, with an initial 
minority of A's, n < Nx* , the fixation of B's is almost 
certain, and we are interested in calculating the exponen- 
tially small probability 0^ = (f)'^{x) that A's fixate. Such 
a probability satisfies the following equation [|[ |[ |j: 

+ r-0li - [r+ + T-^i = o, (i5) 

which is the stationary backward master equation of this 
problem [Q, with boundary conditions = 0?'/'^ = 1- 



At this point, it is convenient to introduce the auxiliary 
quantity 

Vn^cjyi^^-clyi^Vix), (16) 

which is a normalized PDF peaked at x* . From Eq. (|l^) , 
the fixation probability 0^ can be easily obtained, yield- 
ing (Pi = Y.l;}„ Prn- Substituting Eq. (|l^ into Eq. (|l|), 
one arrives at a difference equation for the PDF Vix) 
which reads 

r+ix)P{x)-T-{x)V{x-l/N)=0. (17) 

This equation can be treated with the WKB ansatz 
r{x) = Acce^^^^'''^'^'^'''^- To leading order one 
has T+{x) — T-{x)e^ '^^^ = 0, whose solution is 
S{x) = —S{x) [where S{x) is given by Eq. (||)]. 
In the subleading order, after some algebra, one 
finds Si{x) = (1/2) ln[7+(x)/71(a:)]. Normalizing 
X]n=o^^" — /o 'P(x)dx — 1, and assuming a main 
Gaussian contribution arising from x ~ x*, we find 
^CG - ^/\S" {x*)\/ (2TTN)e-^^^'''\ In the realm of the 
WKB approach, we have thus obtained an expression of 
V{x) that holds for < a; < 1. [It can be checked that 
such a WKB result satisfies Eq. ( p^ ) also near the ab- 
sorbing boundaries From the expression of V{x), 
the fixation probability thus reads 

<t>t = y^pE e^[5(W^)-5(.*)l. (18) 

m=0 

Of special interest is the limit of n <^ Nx* correspond- 
ing to the fixation of a few mutants in a sea of wild-type 
individuals . In this case, it can be shown that Eq. ( [T^ ) 
can be well approximated by (t)'^{x) ~ V{x)/[e^ — 1] 
when w = 0(1), while for small selection A^~^ -C w <C 1, 
(j)^{x) ~ ,jN\S"{x*)\/{2n) Sidy e^^s(y)-S(x')] Q. A 
comparison between theory (Hq) and numerical results, 
using S{x) from Eq. (p^), is shown in Figs. ^ and ^(a) 
and an excellent agreement is observed. 

The WKB theory presented in this section (as well as 
that dedicated to the ACG) is valid as long as w N^^. 
For small selection intensity, w <^ \ the fixation 

probability is often computed using the FPA (or diffu- 
sion approximation) [^|-^, usually considered within the 
linear noise approximation ||^. Thus, for N~^ <^ w ^ 1 
(and > 1), the predictions of the WKB and FPA 
approximations can be compared (together with results 
of numerical simulations) to determine their respective 
domains of validity. 

For the purpose of comparison, it is conve- 
nient to rewrite both WKB and FPA predictions 
in the following form: 0^(a;) ~ ^'(a::)/^(l), 
where *(a;) = Jg" e~ ^o" e(^). From Eq. (|l|), 
one finds that for the WKB approach the expo- 
nent reads 8wkb(^) = A^ In [7+(2;)/7^(z)], whereas 
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FIG. 3: (Color online). The fixation probability 0'* (x) for the 
fMP process: theoretical result (|l^ (solid), numerical calcu- 
lations (dashed) and FPA (dash-dotted), with a = 4,b — 
0.2, c = 0.3, d = 3.8, = 100. Insets: ratio between theoret- 
ical results and those of the FPA, see text, (a) For w = 0.1, 
Nw'^ — 1 and all curves agree well, with an error of about 7% 
in the predictions of the FPA for a; — >■ 0. (b) For w = 0.75, 
Nw^ ^ 1, the curve obtained from the FPA systematically 
deviates from the others and yields exponentially large errors. 
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FIG. 4: (Color online), (a) Fixation probability (^'^(x) as 
function of ty: theoretical result ( [l^ ) (solid), numerical calcu- 
lations (dashed) and FPA (dash-dotted), for o = 1, 6 = 0.2, 
c = 0.3, d — 0.8, and N = 200. (b) Ratio between the predic- 
tions of the FPA and those of our theory vs A'^, for w = 0.25, 
a = 4, 6 = 0.2, c = 0.3, and d = 3.8. The results of the 
FPA deteriorate when both w and A'' increase. In (a) and 
(b), n = 10 thus X = 10/A, and the system follows the fMP. 



for the FPA one has the exponent 0fpa(-z) = 
2Nz [Tiix*) - Tlix*)] I [r+(a;*) + T-{x*)\ §. Hence, 
it can be shown that in the vicinity of x = x* , 0wkb (x) — 
6FPA(a;) ~ Nw^{x - x*)^ |2l|. Therefore, while the 
WKB result ( p^ is accurate for any finite value of w 
[as shown in Fig. ^(a)], the FPA is unable to account 
for fixation and yields exponentially large errors when 
w > N^'^/'^. In fact, the predictions of the FPA (within 
linear noise approximation) are accurate only when the 
selection intensity satisfies w <^ N~-^^^, which is a more 
stringent condition than w <^ 1. This is illustrated in 
Figs. H and ^ which display a comparison between our 
predictions and those of the FPA for various values of w 
and N. 



CONCLUSION 

We have studied fixation in evolutionary games un- 
der non- vanishing selection and elucidated the nontrivial 
relation between selection intensity and effects of demo- 
graphic fiuctuations. This has been achieved by general- 
izing a recent WKB-based theory to account for multi- 
ple absorbing states. This approach naturally accounts 
for non-Gaussian behavior and allows an accurate treat- 
ment of large fluctuations. In the framework of mod- 
els of cooperation dilemmas, we have analytically com- 
puted the QSD (shape of the metastable PDF), MFTs 
and the fixation probabilities beyond the weak selection 
limit. While it does not cover the w — > limit (where 
the FPA holds), our theory agrees excellently with nu- 
merical simulations over a broad range of finite selection 
strength {0 < w < 1), where the FPA generally fails. 
For concreteness, our approach has been illustrated for 
two classes of (formally solvable) 2x2 games, but is nei- 
ther restricted to linear payoffs nor to a specific choice of 
the transition rates . Importantly, our theory can be 
adapted to study evolutionary processes for which there 
is no rigorous analytical treatment (e.g. 3x3 games |^) 
and help generalize the concept of evolutionary stability. 

Acknowledgments: We thank Baruch Meerson for a 
useful discussion. 
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